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Overview 

This  report  describes  the  progress  and  accomplishments  for  the  period  September  25,  2004 
-  June  30,  2008  for  the  MURI  grant  “Accurate  Theoretical  Predictions  of  the  Properties  of 
Energetic  Materials”  (No.  DAAD19-02-1-0176).  This  is  a  multi-university,  comprehensive 
theoretical/computational  research  program  to  develop,  validate,  benchmark,  and  apply  methods 
and  models  that  will  provide  predictive  capabilities  for  energetic  materials.  The  thrust  of  the 
work  is  the  development  of  atomic-level  models  and  ab  initio  quantum  chemistry  methods  that 
are  generally  applicable  to  the  chemical  decomposition  of  condensed-phase  energetic  materials 
under  extreme  conditions. 

The  approaches  include  quantum  mechanics,  molecular  modeling,  Monte  Carlo,  and 
molecular  dynamics  to  yield  state-of-the-art  methods  specifically  designed  for  and  tailored  to 
target  DoD  energetic  materials  research  needs.  The  following  is  a  list  of  the  principal 
investigators,  universities,  and  topics  that  make  up  the  overall  MURI  project: 

•University  of  Missouri  -  Columbia;  Professor  John  E.  Adams  (PI);  Gas-liquid  heat  transfer 
•University  of  Maryland;  Professor  Herman  L.  Ammon  (PI);  Crystal  models 
•University  of  Florida;  Professor  Rodney  J.  Bartlett  (PI);  Quantum  Chemistry 
•North  Carolina  State  University;  Professor  Donald  W.  Brenner  (PI);  Reaction  potentials 
•University  of  Illinois;  Professors  David  Ceperley  and  Richard  Martin  (Pis);  Solid-state 
quantum  mechanics 

•University  of  Missouri  -  Columbia;  Professor  Donald  L.  Thompson  (PI);  Reaction  dynamics 
and  simulations 

•University  of  Minnesota;  Professors  Donald  G.  Truhlar  and  Christopher  J.  Cramer  (Pis); 
Solvation  models 

The  focus  of  the  project  is  on  developing  accurate  methods  for  simulating  physical  and 
chemical  processes  in  condensed-phase  energetic  materials.  The  MURI  team  is  working  closely 
with  researchers  at  DoD  labs,  who  are  contributing  to  the  theoretical  efforts,  providing  data  for 
testing  of  the  models,  or  aiding  in  the  transition  of  the  methods,  models,  and  results  to  DoD 
applications. 

The  major  goals  of  the  project  are: 

•  Models  that  describe  phase  transitions  and  chemical  reactions  in  condensed  materials. 

•  Ab  initio  predictions  of  structures  and  properties  of  solids  at  high  temperatures  and  pressures. 

•  Methods  to  predict  mechanical  properties  and  physical  changes  in  condensed  phases,  including 
thermodynamic  quantities,  phase  diagrams,  vibrational  frequencies,  tensile  strength,  stress/strain 
relations,  and  densities. 

•  Simulation  methods  to  predict  chemical  decomposition  in  condensed  phases,  particularly 
ignition  and  sensitivity  in  response  to  heating  and  shocking. 

•  Methods  for  predicting  temperatures  of  the  condensed  phases  and  flames  resulting  from 
physical  and  chemical  changes,  including  a  predictive  model  for  the  “heat  feedback”  from  the 
flame  to  burning  surface. 

•  Methods  for  predicting  solvation  and  separation  for  energetic  materials  in  supercritical  fluids. 


3 


The  following  are  descriptions  of  the  component  efforts  that  make  up  the  overall  project. 

Structure  and  Solid-State  Density  Prediction 

Herman  L.  Ammon  (Department  of  Chemistry  and  Biochemistry,  University  of  Maryland) 

The  discovery  of  new  energetic  materials  can  be  facilitated,  accelerated  and  made  more 
cost  effective  with  the  use  of  computer  modeling  and  simulations  for  the  identification  of 
compounds  that  have  significant  advantages  over  materials  currently  in  use.  The  quantitative 
estimation  of  key  parameters  can  be  used  to  screen  potential  energetic  candidates  and  permit  the 
selection  of  only  the  most  promising  substances  for  laboratory  synthesis  and  development. 

The  most  significant  properties  or  characteristics  of  a  high  performance  energetic 
material  are  the  molecular  structure  and  elemental  composition,  and  solid-state  heat  of  formation, 
density  and  microstructure  (3-D  structure  with  various  dislocations  and  defects).  Our  primary 
goals  are  the  development  of  procedures  for  the  prediction  of  the  crystal  structures,  accurate 
densities  and  heats  of  formation  of  energetic  materials,  and  to  investigate  the  relationships 
between  crystal  structure/microstructure  and  sensitivity,  compressibility,  polymorphism  and 
crystal  shape. 

A  three  step  model-packing-refmement  procedure  to  predict  the  crystal  structure  of  a 
material  from  an  initial  molecular  model  has  been  developed  and  extended  to  the  major  types  of 
energetics,  bonding  patterns  and  symmetries.  The  MOLPAK  (MOLecular  PAcKing) 
procedure/program  takes  into  account  molecular  structure  and  conformation,  probable  crystal 
packing  arrangements  and  efficiency.  Hypothetical  structures  are  assembled  based  on  the 
structure  of  a  “search  probe”  (molecular  model)  and  molecular  coordination  information 
obtained  from  analyses  of  more  than  a  thousand  experimental  crystal  structures.  The 
coordination  information,  for  the  major  triclinic  to  orthorhombic  space  groups  with  Z’  of  one- 
half  or  one,  covers  approximately  94%  of  the  known  crystal  structures  of  C-H-N-O-F-containing 
materials. 

MOLPAK  utilizes  the  repulsion  part  of  a  6-12  interatomic  potential  to  position  the  search 
probe  elements  in  a  hypothetical  crystal  lattice.  The  repulsion-only  terms  dramatically  accelerate 
the  positioning  calculations  because  the  interatomic  energy  summation  between  a  molecule  and 
its  neighbors  is  terminated  when  the  energy  reaches  a  pre-determined  threshold. 

The  rigid  search  probe  approach  used  in  the  MOLPAK  procedure  to  create  hypothetical 
crystal  structures  is  adequate  if  crystal  forces  do  not  appreciably  alter  the  molecular 
conformation  from  that  in  an  isolated  molecule.  A  procedure  (ROTPAK,  the  ROTational 
PAcKing  program)  has  been  developed  to  allow  conformational  flexibility  concomitant  with 
crystal  packing.  The  addition  of  flexibility  substantially  complicates  the  packing  process  and 
requires  a  fine-tuned  potential  to  balance  the  intra  and  intermolecular  contributions  to  the  overall 
energy.  A  preliminary  example  of  the  potential  of  ROTPAK  is  found  in  packing  calculations  for 
1,3-dinitrocubane.  The  experimental,  MOLPAK-predicted  and  ROTPAK-predicted  densities  are 
1.647,  1.620  and  1.646  g  cm'3,  respectively,  and  the  C-NO2  twist  angles  from  ROTPAK  are 
within  a  few  degrees  of  the  X-ray  values. 

The  third  step  in  the  overall  model-packing-refmement  crystal  structure  determination  is 
refinement  of  the  best  of  the  MOLPAK/ROTPAK  structures  by  optimization  the  unit  cell 
parameters  and  position  and  orientation  of  the  search  probe.  Intramolecular  adjustments  also  are 
necessary  for  ROTPAK  calculations.  The  refinement  entails  minimization  of  the  crystal  lattice 
energy  with  either  the  WMIN/PMIN  or  DMAREL/DMACRYS  lattice  energy  and  refinement 
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procedures/programs.  The  WMIN/PMIN  intermolecular  potential  contains  a  simple  Coulombic 
term  and  standard  6-exp  terms  for  attraction  and  repulsion;  e;  is  a  point  charge  on  atom  i  from 

electrostatic  potential  calculations  (6-3  lg*  basis  set),  rjj  is  the  distance  between  atoms  i  and  j  and 
A,  B  and  C  terms  are  empirical  coefficients.  The  DMAREL/DMACRYS  potential  is  similar, 
but  with  a  distributed  multipole  electrostatic  term  (6-3 lg**  basis  set)  in  place  of  the 
WMIN/PMIN  monopole. 

Ewmin  =  eie/rij "  (AjAj)/rlj  +  (BjBj)exp[-(Cj+Cj)  rjj)] 

The  A-C  coefficients  in  both  potentials  are  obtained  by  fitting  calculated  to  experimental 
structures,  crystal  densities  and  available  heats  of  sublimation.  A  major  effort  was  required  for 
parameterization  of  the  various  atom  functionality.  There  are,  for  example,  different  coefficients 
for  the  nitrogens  in  C-NO2,  N-NO2,  N-NO2  and  ordinary  amines.  The  necessary  program 
infrastructure  for  automatic  coefficient  optimization  (simplex  followed  by  least  squares)  is  in 
place  and  requires  effort  primarily  in  the  set-up  of  input  data  sets  and  experimental  data  quality 
evaluations.  The  limited  experimental  heat  of  sublimation  data  (related  to  lattice  energies) 
possibly  can  be  supplemented  with  ab  initio  energies  with  appropriate  corrections.  Additionally, 
the  solid-state  heats  of  formation  can  be  obtained  by  the  appropriate  combination  of  the 
WMIN/PMIN/DMAREL/DMACRYS  lattice  energies  and  molecular  heats  of  formation  from  ab 
initio  calculations.  Additionally,  these  procedures  and  potentials  can  be  used  to 
investigate/predict  polymorphism  and  crystal  habits. 

Sensitivity,  the  ease  with  which  a  material  undergoes  a  violent  reaction  or  detonation,  is 
defined  in  relation  to  particular  tests  or  stimuli  such  as  impact,  shock,  friction,  thermal  and 
electrostatic  sources  of  energy.  A  number  of  molecular  structure-sensitivity  relationships  have 
been  found  over  the  years  and  good  correlations  observed  within  families  of  energetic  materials. 
Sensitivity  is  perhaps  the  most  complicated  and  least  well  understood  of  the  various  derivative 
properties  of  energetic  materials.  Areas  that  can  be  explored  with  structure  predictions  for 
unknown  materials:  (1)  Density  of  states  relationship  of  Kunz  which  requires  the  crystal 
structure  and  appropriate  ab  initio  crystal  lattice  calculations.  (2)  Possible  relationships  between 
impact  sensitivity  and  bond  strength  plus  lattice  energy.  The  basic  idea  is  that  impact  has  an 
important  frictional  or  heat  component  that  initially  causes  disruption  (melting  or  deformation) 
of  the  crystal  lattice  followed  by  homolytic  cleavage  of  the  weakest  bond  in  the  molecule  leading 
to  detonation.  Lattice  energies  are  available  from  WMIN/PMIN/DMAREL/DMACRYS 
calculations  and  bond  energies  (e.g.,  C-NO2,  N-NO2,  O-NO2)  from  ab  initio  calculations  of  the 
energies  of  the  appropriate  framework  and  nitro  group  radicals.  (3)  A  correlation  between  the 
impact  sensitivity  and  hydrostatic  compressibility  (from  molecular  dynamics  calculations)  for  a 
variety  of  energetics  has  been  observed.  (4)  The  lattice  potentials  developed  can  be  utilized  to 
explore  the  relationship  between  crystal  orientation  and  detonation.  In  PETN,  for  example,  the 
orientation  to  shock  initiation  and  detonation  is  consistent  with  steric  hindrance  to  sheer  at  the 
molecular/crystal  level. 

Details  of  the  MOLPAK,  PMIN,  ROTPAK  and  volume  additivity  codes  follow. 

MOLPAK-1.  The  computer  code  typically  used  to  generate  hypothetical  packing 
arrangements  is  termed  here  MOLPAK-1,  to  distinguish  it  from  a  new  generation  code  under 
development  termed  MOLPAK-2.  MOLPAK  generates  hypothetical  crystal  packing  motifs  that 
mimic  the  experimental  coordination  geometry  information  determined  from  the  analysis  of 
1000’s  of  C-H-N-O-L  crystal  structures  in  the  triclinic,  monoclinic  and  orthorhombic  space 
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groups.  Analysis  of  the  three-dimensional  crystal  packing  arrangements,  with  the 
MOLCON/MOLPAN  (molecular  connectivitv/molecular  packing  analysis)  codes,  revealed  a 
relatively  small  number  of  common  coordination  patterns  and  this  information  was  used  to 
construct  routines  to  handle  54  distinct  patterns.  The  patterns  are  based  on  a  coordination 
number  of  14,  the  number  of  molecules  in  van  der  Waals  contact  with  a  central  molecule  We 
typically  use  only  31  of  the  54  geometries  which  cover  the  majority  of  cases  for  energetic 
materials.  Additionally,  MOLPAK-1  handles  molecules  with  Q,  C2  and  Cs  point  group 
symmetries  in  which  the  molecule  is  coincident  with  appropriate  unit  cell  symmetry.  P-HMX, 
for  example,  has  C,  point  group  symmetry  and  the  molecule  occupies  a  center  of  symmetry  in  the 
P2i/c  space  group. 

In  the  construction  of  hypothetical  crystal  structures  for  each  of  the  coordination  geometries, 
MOLPAK  typically  rotates  the  model  of  interest  (search  probe)  about  three  Eulerian  axes  by 
180°  in  10°  steps  for  a  total  of  6,859  orientations,  each  of  which  is  used  to  generate  a  possible 
crystal  structure.  The  number  is  50,653  for  a  5°  step  search.  MOLPAK  has  been  expanded  to 
save  up  to  60,000  hypothetical  structures  and  to  pass  up  to  5,000  of  the  highest  density  for 
subsequent  lattice  optimization.  Generally,  after  eliminating  identical  solutions  and  averaging 
similar  adjacent  solutions,  the  number  remaining  for  optimization  is  on  the  order  of  1-2,000. 

PMIN.  Each  of  the  MOLPAK-generated  structures  is  refined  by  lattice  energy  minimization 
to  optimize  the  unit  cell  parameters,  orientation  and  position  of  the  search  probe.  The  work¬ 
horse  refinement  code  WMIN,  developed  in  the  late  70 ’s,  has  been  replaced  by  a  simpler,  cleaner 
and  easier  to  understand  code  named  PMIN.  The  major  features  are  (a)  modem  fortran  90/95; 
(b)  step  and  least-squares  optimization;  (c)  facile  multiple  molecule/ fragment  refinement,  such  as 
cation-anion  complexes;  (d)  both  lattice  and  molecular  conformation  refinement;  (e)  ability  to 
easily  use  individual  cross-term  potential  coefficients.  Such  as  for  H-  bonding;  (f)  potential 
energies  can  be  calculated  either  by  Ewald  or  direct  summations  or  a  combination  of  the  two;  (g) 
user  defined  potentials  are  easily  implemented. 

ROTPAK.  Conformational  changes  concomitant  with  the  crystallization  process  can 
occur  in  a  number  of  flexible  molecules.  This  is  evident  if  one  compares  the  structure  of  an 
optimized  model  or  vapor  phase  conformation  (no-neighbor  conformation)  with  that  found  in  the 
solid  state.  The  flexibility  inherent  in  the  eight-membered  HMX  gives  rise  to  crystals  in  which 
the  molecular  point  group  symmetries  are  C,  and  C2.  The  ROTPAK  code,  similar  to  MOLPAK, 
was  developed  to  permit  limited  conformational  changes  to  accompany  the  typical  packing 
process.  Subroutines  have  been  added  to  permit  bond  rotation  (e.g.  C-NO2),  out-of-plane 
bending  (e.g.  C2N-NO2)  and  bond  angle  variations.  The  packing  process  requires  an  appropriate 
mixture  of  the  individual  inter  and  intramolecular  energies.  The  latter  are  currently  evaluated 
with  external  calls  to  Gaussian-03  to  obtain,  for  example,  semi-empirical  energies  such  PM3  or 
ab  initio  b31yp/6-31g*  energies.  An  important  part  of  this  procedure  is  the  lattice  plus 
conformational  refinement  by  total  energy  minimization.  Conformational  refinement  additions 
have  been  added  to  the  step  search  portion  of  PMIN  and  changes  to  the  least-squares  routines  are 
planned. 

The  basic  packing  procedures  in  MOLPAK  and  ROTPAK  differ,  although  both  utilize  the 
same  coordination  geometry  information  for  packing  construction.  MOLPAK  has  a  repulsion- 
only  potential  and  molecules  are  docked  until  a  pre-determined  repulsion  threshold  is  met.  This 
accelerates  packing  because  calculation  of  the  total  repulsion  energy  between  molecules  can  be 
terminated  as  soon  as  the  threshold  is  achieved,  the  distance  is  then  increased  and  the  calculation 
repeated.  In  ROTPAK,  however,  the  intermolecular  energy  between  two  (or  more)  molecules  is 


6 


calculated  with  the  same  q-6-exp  interatomic  potential  used  in  PMIN  and  distances  are  adjusted 
to  a  minimum  energy.  The  results  from  MOLPAK  and  non-flexible  ROTPAK  calculations  are 
similar,  in  most  cases,  and  generally  the  same  after  PMIN  refinement. 

DMAREL  (distributed  multipole  analysis  relaxation).  The  program,  an  alternative  to 
PMIN  for  lattice  refinements,  uses  a  distributed  multipole  electrostatic  term  in  place  of  single 
atom-centered  charges.  The  set-ups  for  the  calculations  are  complicated  and  somewhat  clumsy, 
despite  several  automatic  procedures  we  have  introduced,  and  the  refinements  are  4-5  times 
slower  than  with  PMIN.  Generally,  we  have  preferred  the  easier  and  faster  codes  but  recently 
have  returned  to  the  superior  results  typically  provided  by  DMAREL.  DMAREL  was  developed 
at  University  College  London  and  is  being  replaced  by  DMACRYS. 

VOLUME  ADDITIVITY.  Appropriate  atom  and  group  volumes  can  be  summed  to 
estimate  an  effective  solid-state  molecular  volume  for  a  substance,  which  in  combination  with 
the  molecular  mass  provides  an  estimate  of  the  crystal  density.  With  crystal  structure  data  from 
the  Cambridge  Structural  Database,  four  new  additivity  data  bases  have  been  determined  to 
provide  atom  and  functional  group  volumes  or  densities  for  the  calculation  of  solid  state  densities 
of  single  organic  compounds,  organic  salts  and  multi-fragment  materials  with  the  elements  H,  C, 
B,  N,  O,  F,  P,  S,  Cl,  Br  and  I.  Linear  and  nonlinear  volumes  and  linear  densities  were 
determined  for  108  pre-defined  atoms/functional  groups  from  ~  41,000  crystal  structure  data.  In 
separate  calculations,  the  use  of  atom_codes,  which  do  not  require  a  beforehand  atom  or 
functional  group  specification,  allowed  the  atoms  and  their  connections  to  automatically  define 
the  types  of  atoms  present  in  a  molecule.  This  approach,  applied  to  both  single  molecule  and 
multiple  species  unit  cells  (e.g.,  Z’  >  1,  hydrates,  cation  +  anion),  lead  to  2,805  atomcodes  from 
~  48,000  structures.  With  1,027  crystal  structure  data  not  used  in  the  initial  parameterizations 
and  the  atomcode  volume  parameters,  the  average  percent  difference  between  the  observed  and 
calculated  densities  was  2.3%.  Similarly,  864  new  structures  with  the  108  atom/group  parameters 
gave  a  2.0%  average  volume  difference  with  the  linear  volume  parameters.  This  work  has  been 
published  in  Propellants,  Explosives,  Pyrotechnics. 

MOLPAK-2.  The  present  MOLPAK- 1  code  is  based  on  constructing  hypothetical  crystal 
packing  arrangements  to  mimic  experimental  crystal  structure  motifs.  MOLPAK- 1  has  certain 
limitations,  not  the  least  of  which  is  the  inability  to  handle  new  crystal  symmetries  and  space 
groups  that  previously  have  not  been  coded.  MOLPAK-2  is  designed  to  determine  the  most 
likely  packing  arrangement  from  little  more  than  the  structure  of  the  search  probe  and  crystal 
symmetry  information,  such  as  the  space  group  and  symmetry  operations.  The  starting  point  for 
MOLPAK-2  involves  finding  low  energy  pairs  of  molecules  that  are  related  by  unit  cell 
translations,  inversion  centers,  two-fold  screw  axes,  mirror  glide  planes  and  other  possible 
symmetries  dictated  by  the  unit  cell  symmetry.  These  motifs  of  molecules  are  expanded  to 
possible  unit  cells.  Unlikely  solutions  are  culled  by  the  elimination  of  close  intermolecular 
contact  structures  and  accepting  only  those  with  densities  within  0.1  g/cc  of  a  pre-estimated, 
volume  additivity  values.  Complete  lattice  energy  calculations  with  PMIN  are  used  to  select  the 
most  probable  solutions. 

The  codes  for  the  triclinic  space  groups  PI  and  P-1  and  the  monoclinic  P2i  are  complete 
and  have  been  tested.  Work  is  progressing  on  the  considerably  more  complicated  (and  potentially 
important)  P2i/c.  A  possible  enhancement  of  the  basic  MOLPAK-2  concept  is  the  simultaneous 
use  of  information  from  two  types  of  pairs. 
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First  Principles  Simulations  of  Energetic  Materials 

David  Ceperley  and  Richard  Martin  (Department  of  Physics,  University  of  Illinois,  Urbana- 
Champaign) 

A  large  part  of  our  effort  has  gone  into  developing  the  quantum  Monte  Carlo  code 
QMCPACK  to  allow  more  efficient  QMC  calculations  of  energetic  materials.  The  code  has  been 
publicly  released.  A  presentation  of  Monte  Carlo  methods  and  the  use  of  the  code  was  made  at 
the  Army  Research  Lab  in  July,  2007. 

In  our  developments  of  new  methods,  we  extended  QMC  methods  to  compute  energy 
differences.  This  is  an  important  general  problem  in  QMC  and  we  are  developing  methods  for  an 
efficient  and  accurate  determination  of  dissociation  energies,  reaction  barriers  and  the  entire 
energy  curve  for  the  bond-dissociation,  and  weak  intermolecular  interactions  in  molecular  crystal 
and  van  der  Waals  complexes. 

These  problems  are  especially  important  in  energetic  materials  and  QMC  will  allow  us 
to  compute  such  properties  for  as  big  complexes  as  the  benzene  and  RDX  dimers.  Our  approach 
based  upon  the  reptation  Monte  Carlo  and  is  being  tested  on  dissociation  of  diatomic  molecules. 

We  have  completed  work  on  a  new  method  to  determine  the  finite  size  error  on  the 
energy  for  quantum  systems  in  periodic  boundary  conditions.  This  is  important  for  work  on 
condensed  phases  and  we  have  tested  the  method  on  the  uniform  electron  gas,  solid  silicon  and 
solid  nitrogen.  The  results  of  this  work  are  published  in  Phys.  Rev.  Letters. 

We  are  developing  coupled  electronic-ionic  Monte  Carlo  methods  for  applications 
beyond  hydrogen.  Preliminary  testing  for  systems  of  water  molecules  is  in  progress.  This  is  a 
long  range  project  that  has  the  potential  to  provide  a  new  level  of  accuracy  and  reliability  of 
predictions  for  materials  under  extreme  conditions. 

We  have  initiated  work  to  implement  and  apply  recently-developed  optimization 
methods  in  Variational  Monte  Carlo  (VMC)  which  the  accuracy  and  allow  atomic  geometries 
can  be  optimized  simultaneously  with  minimizing  the  parameters  in  the  electronic  wave 
function.  The  main  goal  is  to  treat  using  methods  that  have  recently  been  applied  to  benzene  and 
benzene  pairs.  The  first  application  is  to  benchmark  calculations  for  a  system  like  H2  adsorbed 
on  benzene;  a  difficult  case  with  very  binding  energy. 


Ab  initio  Predictions  for  Potential  Energy  Surfaces  for  Chemical  Reactions 
Andrew  G.  Taube,  Josh  McClellan,  and  Rodney  J.  Bartlett  (University  of  Florida,  Quantum 
Theory  Project,  Gainesville,  FL  32611-8435) 

Decomposition  Pathways  of  Nitroethane 

The  2006  report  summarized  this  extensive  work,  so  that  will  not  be  repeated  here.  Suffice 
it  to  say,  that  the  major  new  development  was  the  introduction  of  frozen  natural  orbitals  (NO) 
into  coupled-cluster  (CC)  theory,  which  have  the  effect  of  reducing  the  computational 
dependence  of  the  CC  calculation  by  about  an  order  of  magnitude  with  no  loss  in  accuracy.  This 
required  that  entirely  new  theory  for  the  analytical  derivatives  (forces)  of  NO  based  CC  be 
developed  and  implemented  into  ACES  II.  Using  this  new  capability,  a  very  thorough,  predictive 
study  of  the  decomposition  and  rearrangements  of  the  prototype  energetic  nitroethane  molecule 
was  presented  with  state-of-the-art  CC  methods.  The  discrepancies  with  the  results  from  DFT 
calculations  were  considerable,  suggesting  that  DFT,  as  used  for  these  systems,  is  not  a  reliable 
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tool  to  understand  the  comparative  importance  of  the  various  reaction  paths.  The  lowest  path  in 
the  case  of  nitroethane  was  found  to  be  the  HONO  elimination.  This  paper  is  now  published  in 
the  J.  Chemical  Physics. 

ACCSD(T):  The  new  “gold-standard’? 

For  some  years,  the  gold  standard  for  quantum  chemistry  has  been  the  CCSD(T)  method, 
as  it  is  fairly  widely  applicable  to  moderate  sized  molecules.  However,  the  (T)  part  means  that 
coupled-cluster  single  and  double  excitations  (CCSD)  is  augmented  by  a  non-iterative, 
perturbative,  triples  correction.  It  is  well-known  in  the  field,  that  once  bonds  start  to  be  broken, 
as  in  a  transition  state  (TS),  that  such  perturbative  corrections  can  fail  dramatically,  particularly 
when  using  a  spin-adapted  RHF  reference  function.  An  illustration  is  shown  below  for  breaking 
the  LiF  bond. 


Lithium  Fluoride 


F:  (9s  6p  Id]  ->  [4s  3p  Id]  basis,  Is, 2s  frozen 
From  Bauschlicher,  Langhoff  JCP  89,  4246  (1988) 

DMRG  from  Legeza,  Roder,  Hess,  Mol.  Phys.  101,  2019  (2003) 


What  should  be  observed  from  this  figure  is  that  at  larger  internuclear  separation,  that 
CCSD(T)  falls  precipitously  off  the  correct  curve,  here  defined  by  the  full  Cl.  Note,  however, 
that  once  the  new  ACCSD(T)  method  is  applied  to  this  problem,  that  a  good  curve  is  obtained  all 
the  way  to  the  dissociation  limit. 

This  is  not  an  isolated  situation,  as  this  can  be  observed  in  nearly  all  bond  breaking  as 
illustrated  for  the  triply  bonded  CO  molecule. 
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The  failure  of  CCSD(T)  as  the  bond  length  is  extended  beyond  its  equilibrium  geometry  suggests 
that  for  transition  states  that  are  farther  from  the  molecule’s  equilibrium  geometry,  it  will  be 
difficult  to  be  definitive  for  TS  geometries  and  their  associated  activation  barriers.  However,  as 
we  have  already  seen  in  the  nitroethane  example,  DFT  cannot  be  expected  to  sort  out  competitive 
paths  correctly,  so  there  is  a  need  for  the  best  CC  methods  to  be  applicable  to  the  problem. 

The  ACCSD(T)  method  is  still  pertubative,  like  CCSD(T),  which  gives  it  its  high 
applicability;  but  unlike  that  method,  benefits  from  the  use  of  a  left-hand  eigenvector,  A,  that  is 
determined  prior  to  imposing  the  (T)  correction.  By  so  doing,  many  of  the  effects  of  the 
perturbation  approximation  are  removed,  providing  a  much  more  stable  treatment  of  molecular 
geometries  away  from  equilibrium.  A  case  in  point  is  the  ‘triple-whammy’  fragmentation  of 
RDX  to  three  equivalent  units.  This  has  been  supported  by  Yuan  Lee’s  molecular  beam 
experiments,  but  is  disputed  by  most  explosives  experts  who  anticipate  that  there  is  a  direct  bond 
fission  of  the  nitramine  bond.  Some  of  the  issues  include  the  comparative  relevance  of  the 
unimolecular  dissociation  path  versus  that  observed  in  condensed  phases,  where  bimolecular 
mechanisms  might  predominate;  and  the  comparative  entroies  along  the  different  unimolecular 
paths.  This  feature  points  to  a  need  to  obtain  the  more  relevant  free  energy  paths,  at  the  cost  of 
doing  a  vibrational  analysis  for  the  zero-point  and  entropy  at  every  relevant  geometry  along  the 
dissociation  path.  We  are  currently  trying  to  help  resolve  this  question  by  using  the  ACCSD(T) 
method  for  those  TS’s. 

As  in  all  ab  initio  studies  that  require  searching  potential  energy  surfaces  to  locate  TS’s 
and  equilibrium  geometries,  the  development  of  analytically  computed  gradients  is  requisite.  So 
much  of  our  time  was  necessarily  devoted  to  this  development,  that  is  particularly  difficult  for 
ACCSD(T).  However,  this  has  now  been  done  successfully  and  we  have  applied  the  theory  to  a 
variety  of  TS  studies.  This  work  has  been  published  in  two  papers  in  the  J.  Chemical  Physics. 
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Adaptive  Ab  Initio 

The  basic  reason  for  developing  the  adapted  ab  initio  theory  (AAT),  is  to  offer  an 
alternative  to  the  current  semi-empirical  quantum  chemical  approaches  for  very  large  scale, 
quantum  mechanically  based  forces  for  classical  molecular  dynamics  (MD)  simulations  of 
condensed  phase  explosives.  The  more  sophisticated  QM  methods,  like  CC  theory,  MBPT,  and 
even  DFT,  cannot  be  applied  on  a  rapid  enough  time  scale  to  tie  naturally  with  classical 
molecular  dynamics  methods.  Yet  we  require  this  capability  to  allow  for  bonds  to  be  broken  in 
detonation  while  we  use  MD  to  observe  the  time  development  of  the  species.  The  alternative 
today,  is  to  employ  classical  force  fields  for  the  MD  step,  but  such  approaches  have  to 
compromise  on  the  details  of  the  reactions  and  when  applicable,  on  the  spectroscopic  signatures 
that  might  be  observed  for  the  phenomena.  Hence,  only  semi-empirical  theory  like  MNDO, 
AMI,  PM3,  and  their  many  variants  can  be  used  today.  Yet  all  these  are  easily  shown  to  be 
unable  to  correctly  describe  the  potential  energy  surfaces  required  for  bond  breaking. 

The  philosophy  of  AAT  is  to  use  rigorous,  ab  initio  theory  based  upon  a  combination  of 
wavefunction  and  density  functional  theory  as  the  foundation  for  a  new,  very  efficient  QM 
evaluation  of  energies  and  structures  for  -1000  atoms.  All  approximations  introduced — of  which 
none  are  semi-empirical — are  assessed  compared  to  accurate  correlated  theory,  to  determine  their 
viability.  Yet  at  the  same  time,  these  approximations  are  made  with  the  objective  of  routinely 
treating  -1000  atoms.  Currently,  this  approach  is  being  put  into  ACES  III,  the  new  parallel 
version  of  ACES,  to  make  very  large-scale  applications  and  to  further  assess  the  strengths  and 
weaknesses  of  AAT,  to  identify  how  it  can  be  improved. 

Plans  For  The  Future 

We  have  now  shown  that  ACES  III  scales  at  100%  efficiency  up  to  512  processors  (see 
below),  which  is  all  we  have  been  able  to  access  at  the  DoD  HPCMO  installations.  In  the 
transition  to  petascale  applications,  of  which  DoD  is  a  major  player,  this  has  to  be  fixed  for 
future  developments.  Only  with  routine  access  to  >1000  processors  for  relatively  short  periods  of 
time  for  debugging  purposes,  like  15  minutes,  can  the  software  be  adequately  developed  to  begin 
to  explore  even  the  terascale  on  the  way  to  the  petascale. 

With  the  development  of  ACES  III,  which  has  been  written  from  scratch  even  to  writing 
an  entirely  new  integral  and  integral  derivative  program;  we  are  poised  to  make  some  of  the 
definitive  computations  of  great  interest  to  ARO.  An  example  is  the  RDX  molecule’s 
unimolecular  dissociation  paths,  where  with  ACEs  Ill’s  analytical  gradient  methods,  this  path 
can  be  studied  in  depth.  Even  the  paths  for  the  dimer  of  RDX  as  it  falls  apart  can  be  described. 
For  rigid  RDX  molecules,  MANY  points  have  been  done  by  Szalewicz  and  Rice  to  fit  a  force 
field,  which  is  already  a  monstrous  calculation,  but  the  real  dissociation  paths  demand  relaxing 
its  geometry  in  all  possible  ways.  Such  a  study  has  been  well  beyond  the  capability  of  nay  prior 
CC  work,  but  can  now  be  envisioned  with  ACES  III  and  adequate  computer  resources. 

Furthermore,  ACES  III  also  offers  excited  states  with  the  equation-of-motion  (EOM- 
CC)  method,  another  unique  feature  in  such  a  parallel  CC  program  (see  below  for  scaling). 
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Luciferin  CCSD(T) 


Jul  17, 
08 


•  CCSD  on  128  processors 

•  One  iteration:  23  min 

•  Total  12  iterations:  275  min 

•  (T) 

•  Hardest  8  occupied  orbitals:  420  min  on  128  processors 

•  Total  48  correlated  orbitals:  420  min  on  768  processors 
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Succrose  CCSD  scaling 
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Succrose  CCSD  super  linear  scaling 


Jul  17, 
08 


•  CCSD  iteration 

•  32  processors  909  min 

•  512  processors  24  min,  ideal:  57  min 


36 


ACES  III 


Time  for  one  CCSD-EOM  iteration 


4.4  minutes 
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Progress  toward  the  Development  of  Transferable  Bond  Order  Potentials  for 
CHNO 

Donald  W.  Brenner  (North  Carolina  State  University) 

We  have  developed  analytic  potential  functions  for  CHNO  clusters,  solids,  and  molecular 
materials  that  have  a  transferable  functional  form  that  is  independent  of  molecular  structure.  The 
bonding  formalism  is  based  on  a  second-moment  approximation  in  which  a  bond-order  whose 
value  is  proportional  to  the  inverse  square  root  of  local  coordination  is  coupled  to  an  attractive 
pair  term  that  represents  bonding  from  valence  electrons.  Together  with  a  pair-additive,  largely 
repulsive  term,  this  functional  form  is  capable  of  reproducing  Pauling  bond  length/bond 
order/bond  energy  relations  over  a  wide  range  of  coordinations,  while  allowing  bond  breaking 
and  forming  with  appropriate  effective  atomic  rehybridation  in  a  computationally  efficient 
manner.  These  functions  represent  an  important  bridge  between  first  principles  methods  and 
large-scale  atomic  simulations. 

Determination  of  functional  forms  and  parameter  fitting  for  transferable  pair  terms,  bond 
order  and  screening  functions,  as  well  as  a  novel  molecular  topology  function  were  all  completed 
and  incorporated  into  molecular  dynamics  codes  at  both  NC  State  and  the  ARL.  Part  of  the 
fitting  process  for  oxygen  made  extensive  use  of  first  principles  results  generated  within  the 
MURI  team.  Extensive  testing  of  the  potential,  including  comparison  to  additional  first  principles 
results,  has  revealed  that  it  is  robust  with  respect  to  reactive  dynamics  and  the  production  of 
intermediate  compounds  that  were  not  an  explicit  part  of  the  fitting  data  base. 

While  the  development  of  the  functional  form  and  the  delivery  of  a  working  code  to  the 
ARL  was  successful,  testing  has  revealed  the  functional  form  and  its  implementation  to  be  too 
slow  for  many  of  the  uses  targeted  by  ARL  researchers.  Over  the  last  year  we  have  worked  to 
make  the  potential  and  associated  code  more  computationally  efficient  without  loss  of  accuracy 
or  transferability.  As  an  example,  the  screening  function  in  the  topology  expression  was  replaced 
with  a  distance-dependent  function,  and  a  new  implementation  of  the  screening  function  that 
takes  advantage  of  inherent  properties  of  the  function  was  accomplished.  More  efficient  neighbor 
indexing  was  also  implemented  within  our  codes.  These  changes  have  sped  the  code  up  over  a 
factor  of  60,  which  broadens  its  use  for  DoD  and  related  applications. 


Flame-Surface  Energy  Feedback 

John  E.  Adams  (University  of  Missouri  -  Columbia) 

Over  the  period  of  this  project,  we  completed  the  characterization  of  the  basic  features  of 
energy  feedback  from  hot  combustion  products  to  the  liquid  layer  atop  a  deflagrating  propellant 
in  simulations  of  both  prototypical  model  systems  and  actual  energetic  species.  In  all  the 
systems  investigated,  we  have  found  that  that  the  net  transfer  of  energy  to  the  liquid  surface  is 
captured  in  a  simple  kinematic  model,  the  general  form  of  which  was  devised  originally  by  Baule 
to  describe  gas  scattering  from  a  solid  surface.  This  model  previously  was  identified  by 
Nathanson  as  providing  a  useful  depiction  of  his  experimental  measurements  of  energy  transfer 
at  a  liquid  interface,  but  we  have  been  able  to  test  the  model  over  a  range  of  conditions  wider 
than  that  accessible  in  the  experiments  and  to  verify  in  more  detail  its  general  utility  in 
describing  the  net  flow  of  energy  back  to  the  liquid  surface.  We  also  have  assessed  the  role 
played  by  the  internal  vibrational  degrees  of  freedom  of  a  molecular  liquid  in  collisional  energy 
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transfer  and  have  determined  the  extent  to  which  the  initial  translational  energy  of  an  impinging 
gas  molecule  is  deposited  in  rotation  of  that  species.  The  investigations  yielding  these 
generalizations  are  described  in  more  detail  below. 

Our  simulation  methodology  was  calibrated  in  studies  of  energy  exchange  in  a  simple, 
prototypical  atomic  system  (Lennard-Jones  interactions).  Having  first  equilibrated  a  sample  of 
the  neat  liquid  under  conditions  of  constant  density  and  temperature,  we  expanded  the  simulation 
cell  along  one  dimension  thereby  creating  two  interfaces.  (The  potential  energy  function  cutoffs 
are  such  that  the  two  interfaces  do  not  “see”  one  another.)  The  resulting  system  is  again 
equilibrated,  and  a  series  of  temporally  well-spaced  system  snapshots  is  collected  for  use  in 
setting  the  initial  conditions  of  the  subsequent  microcanonical  trajectories  from  which  final 
energy  spectra  and  trapping  probabilities  are  extracted. 

For  this  system,  the  significant  trapping  probabilities  obtained  attest  to  the  high 
efficiency  of  energy  transfer  for  a  wide  range  of  projectile/target  combinations.  The  figure  here 
corresponds  to  results  generated  using  a  high  initial  incident  energy,  but  even  here  we  find  that  a 
sizeable  fraction  of  the  incident  species  lose 
sufficient  energy  (on  average  the  energy  loss  is  over 
70%;  the  thermal  accommodation  coefficient  is, 
roughly,  0.8)  that  they  are  (essentially)  equilibrated 
with  the  surface.  Of  particular  importance  in  these 
results,  though,  is  the  dependence  on  ju,  the  ratio  of 
the  incident  particle  mass  to  the  surface  atom  mass, 
with  the  maximum  in  the  trapping  probability 
appearing  when  the  particles  are  equal  in  mass  and 
with  a  suppression  of  the  sensitivity  of  the  trapping 
to  n  at  higher  surface  temperatures  (the  upper  set  of  points). 

One  also  expects  that  the  probability  of  becoming  trapped  to  depend  on  the  strength  of 

the  attractive  interaction  between  the  incident  and  surface 
species.  Indeed,  that  expectation  was  borne  out  in  the 
calculations  summarized  in  the  next  figure  in  which  we  varied 
the  interaction  well  depth  (either  by  scaling  the  full 
interaction  energy — the  solid  lines — or  by  modifying  only  the 
attractive  component  of  that  energy,  thereby  retaining  the  full 
short-range  repulsive  component — the  dashed  lines).  It  is  not 
surprising  here  that  a  deeper  attractive  well  leads  to  increased 
trapping  or  that  more  energy  is  transferred  (and  thus  more 
trapping  occurs)  when  the  full  short-range  repulsive 
interaction  is  included  in  the  potential,  but  it  is  again 
interesting  to  note  the  strong  surface  temperature  dependence 
of  the  trapping  probability.  Trapping  on  the  cooler  surface, 
while  less  probable  for  small  attractive  well  depths,  actually  becomes  more  probable  than 
trapping  at  the  warmer  surface  when  the  well  depth  is  relatively  large.  Again,  we  find  a  notable 
suppression  of  the  sensitivity  of  the  trapping  (and  hence  of  the  energy  transfer)  to  the  initial 
collision  conditions  when  the  surface  temperature  is  increased. 

The  real  contribution  to  our  understanding  of  energy  feedback  in  flames  comes, 
however,  from  recognition  that  these  results  evince  a  quite  simple  model  of  the  underlying 
particle  dynamics,  one  that  derives  from  the  work  of  Baule  dealing  with  gas-solid  surface  energy 
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transfer.  If  the  incident  atom  is  assumed  to  be  normally  incident  with  kinetic  energy  Et  on  the 
final  atom  in  a  one-dimensional  harmonic  atomic  chain,  then  one  can  show  that  the  fractional 
energy  transfer  is  given  by  the  expression 


A E 
E, 


4// 


|l-  cosj(l 


2  •  2 

/u  sin 


\V2  .  2 

X)  + / u  sin  x 


■][i+(r-2*r,  )/£,], 


where  //  again  is  the  incident/target  atomic  mass  ratio,  %  is  a  deflection  angle  (which  enters 
because  the  line-of-centers  of  the  initial  collision  may  not  coincide  with  the  surface  normal),  V  is 
the  depth  of  the  gas-surface  interaction  potential  (assumed  to  be  a  square  well),  and  Tt  is  the 
temperature  of  the  bulk  liquid.  This  same  expression  was  used  by  Nathanson  in  interpreting  the 
general  trends  in  his  experimental  studies,  but  by  extending  the  range  of  system  parameters 
studied  beyond  those  currently  tractable  by  experiment,  we  can  show  that  this  estimate  of  the 
energy  transfer  has  broad  applicability  to  energy  transfer  at  liquid  surfaces.  It  indeed  predicts 
that  energy  transfer  will  increase  as  ju  approaches  a  value  of  one  and  that  the  maximum  energy 
transfer  will  occur  at  that  point,  although  in  fairness  it  should  be  noted  that  the  dynamical 
assumptions  going  into  the  derivation  of  the  kinematic  equation  begin  to  fail  when  the  mass  ratio 
exceeds  unity.  We  also  have  shown  that  it  is  capable  of  producing  results  that  agree 
quantitatively  with  those  deriving  from  simulation. 

Where  the  kinematic  model  devised  for  gas-solid  inelastic  scattering  misses  the  mark  in 
describing  gas-liquid  energy  transfer  is  in  its  prediction  of  the  surface  temperature  dependence. 
For  high  incident  energies  (the  case  relevant  to  the 
data  depicted  thus  far),  the  surface  temperature  004 

correction  appearing  in  the  final  factor  of  the  " 135 

kinematic  equation  becomes  negligible,  and  yet  we  o  * 

find  a  distinctly  non-negligible  dependence  of  002 

trapping  on  the  liquid  temperature.  An  examination  0.015 

of  the  average  local  density  of  the  simulated  lamella  o.oi 

as  a  function  of  distance  along  the  axis  of  the  0005 

simulation  cell  that  lies  perpendicular  to  the  average  0 

surface  planes  makes  it  clear  that  the  liquid  surface  z  <A> 

roughens  substantially  with  increasing  temperature.  Not  only  does  the  width  of  the  surface 
region  increase,  but  the  position  of  the  Gibbs  surface  advances  farther  from  the  center  of  the 
lamella.  Thus,  at  the  higher  temperatures,  an  incident  particle  is  more  likely  to  be  scattered 
inelastically  away  from  the  specular  direction,  more  likely  to  experience  multiple  collisions  with 
surface  species,  and  more  likely  (due  to  an  increase  in  the  vapor  pressure  of  the  liquid)  to 
encounter  other  gas  atoms  in  the  near-surface  region,  either  before  or  after  collision  with  surface 
atoms.  With  an  increase  in  the  average  number  of  collisions  between  the  incident  and  surface 
atoms,  one  anticipates  that  the  incident  atoms  will  lose  their  memory  of  the  initial  collision 
conditions  and  that  trapping  therefore  will  be  less  sensitive  to  factors  such  as  the  depth  of  the 
interaction  well  or  the  mass  ratio  of  the  collision  partners.  This  expected  behavior  is  precisely 
what  we  see  in  our  simulations. 

Having  established  that  for  a  simple  model  system,  kinematics  and  temperature- 
dependent  surface  roughening  largely  govern  the  observed  energy  transfer,  we  shifted  our  focus 
to  molecular  liquids,  specifically  to  nitromethane.  Initialization  of  the  system  was  performed  in 
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the  same  way  as  for  the  Lennard-Jones  atomic  liquid,  with  the  assumed  potential  energy  function 
for  nitromethane  being  that  described  by  Sorescu,  Rice,  and  Thompson. 


melt,  expand  cell 
in  z-direction  (3-d 
periodic  boundary 
conditions) 

- ► 


From  the  average  number  density  of  nitromethane  molecules  in  our  simulation  cell  following  the 
creation  of  the  free  surfaces  and  subsequent  equilibration,  we  can  easily  determine  the 
equilibrium  vapor  pressure  above  the  liquid.  The  value  obtained  at  360  K,  which  is  0.62  bar, 
agrees  well  with  the  value  (0.66  bar)  estimated  from  literature  values  for  the  enthalpy  of 
vaporization.  Although  the  potential  function  adopted  for  nitromethane  was  optimized  to 
describe  the  crystalline  phase  accurately,  that  function  clearly  is  providing  a  very  reasonable 
description  of  our  liquid-vapor  equilibrium  as  well. 

We  take  the  impinging  atoms  (with 

masses  taken  to  be  that  of  argon)  to  be  incident  Ar  Scattering  from  Nitromethane:  12.5  kJ/mol 

'  incident  energy 

with  a  kinetic  energy  of  4.2RTe.  With  this  initial 
energy,  35%  of  the  atoms  become  trapped  at  the 
surface  of  the  liquid — the  trapping  fraction  here 
is  essentially  the  same  as  that  obtained  in  the 
Lennard-Jones  case  for  the  equivalent  mass  ratio 
and  surface  roughness — but  even  those  that  are 
not  trapped  are  well  on  their  way  to  becoming 
thermalized  with  the  surface.  The  distribution  of 
scattered  atom  energies  shown  in  the  figure  to 
the  right  translates  into  a  thermal  accommodation  coefficient  of  0.81,  which  is  very  similar  in 
magnitude  to  the  values  calculated  in  the  case  of  the  simple  atomic  liquid  described  previously. 

It  has  also  proven  useful  to  examine  a  simplified  model  of  nitromethane,  consisting  of 
structureless  methyl  and  nitro  groups  bound  by  the  same  Morse  potential  used  to  describe  the 
carbon-nitrogen  bond  in  the  full  nitromethane  molecule.  Our  initial  interest  in  this  diatomic 
model  stemmed  from  wondering  whether  we  could  collisionally  induce  a  decomposition  reaction 
at  a  liquid  surface.  (Given  the  paucity  of  full  reactive  potential  functions  that  describe 
dissociation  to  the  correct  structurally  relaxed  products,  we  believe  that  the  diatomic  model 
system  at  least  provides  a  test  of  whether  dissociation  might  be  feasible  under  realistic 
deflagration  conditions.)  However,  even  for  initial  incident  energies  well  in  excess  of  the  C-N 
bond  energy,  no  evidence  of  dissociation  was  observed.  This  result  suggests  that  energy 
transferred  to  translation  or  rotation  of  the  surface  species  is  rapidly  lost  to  the  bulk  and  that  the 
probability  that  the  surface  species  will  undergo  a  collision  that  leads  to  direct  excitation  of  the 
diatomic  above  its  bond  dissociation  threshold  is  very  low.  Nonetheless,  energy  transfer  from 


■  All  incident  atoms 
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the  incident  atoms  to  the  surface  in  this  case  is  quite  efficient.  Whether  we  take  the  attraction 
between  the  incident  atom  and  the  diatomic  to  be  strong  or  absent  entirely,  we  find  that  the 
thermal  accommodation  coefficient  is  greater  than  0.60.  Certainly,  the  fraction  of  the  initially 
incident  atoms  being  trapped  by  the  surface  changes  as  that  attraction  is  varied,  but  in  all  cases 
the  majority  of  the  incident  kinetic  energy  is  deposited  in  the  surface. 

To  understand  our  results  for  the  full  nitromethane  simulation,  we  again  fall  back  on  the 
kinematic  model  that  described  the  Lennard-Jones  system  so  well.  It  is  not  quite  as  easy  to  know 
how  to  evaluate  some  of  the  constants  appearing  in  the  equation  for  the  fractional  energy  transfer 
in  this  case.  The  second  factor  in  the  equation,  which  accounts  for  collisions  that  are  “off- 
center”  clearly  is  problematic  when  the  target  is  a  molecule,  and  the  final  factor,  which  includes 
the  attractive  well  depth,  does  not  account  for  the  differences  in  the  interactions  when  various 
sites  on  the  nitromethane  are  hit.  If,  however,  we  estimate  the  angular  factor  to  be  approximately 
unity  on  average  for  normal  atom  incidence  and  if  we  assume  that  the  attractive  well  depth  is 
negligible,  then  the  equation  predicts  a  fractional  energy  transfer  at  360  K  of  0.50,  which  agrees 
very  well  with  the  value  of  0.52  that  we  obtain  from  a  final  energy  analysis  of  the  scattered 
atoms.  If  we  compare  the  results  obtained  for  an  incident  energy  about  seven  times  higher, 
corresponding  to  the  incident  energy  used  in  the  Lennard-Jones  system  simulations,  then  both  the 
kinematic  equation  and  our  simulations  predict  a  fractional  energy  transfer  of  0.89.  These  recent 
results  thus  suggest  that  the  most  straightforward  way  of  grafting  detailed  molecular  information 
onto  extant  continuum  combustion  models  would  be  to  focus  just  on  the  kinematics,  that  the 
mass  ratios  and  simple  energetic  corrections  to  the  collision  events  capture  the  essentials  of  the 
dynamics. 


For  larger  molecules  of  interest  as 
energetic  materials,  we  recognize  that  we 
need  to  consider  the  possibility  of  energy 
transfer  to  the  internal  vibrational  degrees  of 
freedom  of  these  molecules  by  direct 
collision  rather  than  by  general  heating  of  the 
liquid.  Given  the  role  of  kinematics  in  the 
energy  transfer,  however,  we  speculated  that 
the  vibrational  degrees  of  freedom  of  the 
nitromethane  molecules  play  at  most  a  minor  role  in  the  collisional  energy  transfer.  To 
investigate  this  issue,  we  conducted  a  series  of  simulations  in  which  the  nitromethane  molecules 
were  “frozen”  at  their  equilibrium  geometries. 


Z(A) 


(The  SHAKE  algorithm  was  applied  in  order  to 
constrain  the  molecular  motion.)  The  density 
profile  of  the  lamella  stemming  from  this 
constrained  system  is  shown  in  the  plot  below 
along  with  the  corresponding  profile  of  the 
unconstrained,  “flexible”  liquid.  We  found 
very  little  statistical  difference  between  the  two 
profiles,  certainly  not  enough  to  invalidate  a 
comparison  of  the  energy  transfer  in  the  two 
systems.  The  distribution  of  final  argon  atom 

energies  in  the  rigid-nitromethane  case  is  shown  below.  Qualitatively,  there  is  very  little 
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difference  between  the  distribution  of  scattered  atom  energies  shown  here  and  the  distribution 
obtained  using  the  fully  flexible  nitromethane  force  field.  We  find  a  somewhat  greater 
propensity  for  trapping  in  the  rigid  system  that  perhaps  stems  from  the  absence  of  a  channel  for 
feedback  of  energy  from  excited  vibrational  modes  of  the  surface  molecules  to  the  argon  atoms. 
Nonetheless,  energy  transfer  remains  efficient  in  this  system.  The  estimated  thermal 
accommodation  coefficient  in  the  rigid-nitromethane  case  is  0.89,  which  is  within  10%  of  the 
value  obtained  when  using  the  flexible-nitromethane  model.  That  the  surface  appears  to  be 
somewhat  “stiffer”  when  the  nitromethane  molecules  are  held  rigid  is  also  consistent  with  a 
small  preference  for  forward  scattering  of  the  argon  atoms,  whereas  the  distribution  of  scattering 
angles  is  essentially  symmetric  in  the  forward  and  backwards  directions  when  the  rigidity 
criterion  is  relaxed.  Overall,  however,  we  conclude  that  internal  vibrational  degrees  of  freedom 
play  a  negligible  role  in  the  net  gas-surface  energy  transfer  at  a  liquid  nitromethane  surface. 

Finally,  we  have  extended  our  studies  of  collisional  energy  transfer  between  hot  gas 
atoms  and  atomic  and  molecular  liquid  surfaces  by  considering  the  collisions  of  simple 
molecular  species  with  a  prototypical  atomic  liquid  surface.  In  our  previous  work  the  role  that 
rotational  excitation  of  the  impinging  species  might  play  in  a  general  gas-surface  collision  event 
was  left  undetermined.  To  examine  this  last  process,  we  are  simulating  collisions  of  rigid  CO2 
molecules  with  a  Lennard-Jones  atomic  liquid.  These  incident  molecules,  while  translationally 
“hot”,  are  taken  to  be  nonrotating  initially.  Thus,  any  rotational  energy  present  in  the  reflected 
species  thus  is  necessarily  attributable  to  the  collision  event. 

In  the  upper  panel  of  the  figure  shown  below,  we  display  distributions  of  the  final 
rotational  energies  (given  as  a  fraction  of  the  initial  kinetic  energy)  of  scattered  CO2  molecules, 
each  corresponding  to  a  different  incident-molecule  kinetic  energy.  Note  that  each  of  these 
distributions  indicate  that  a  substantial  fraction  of  the  incident  energy  ultimately  is  deposited  in 
the  rotational  degrees  of  freedom  of  the  scattered  species,  with  this  energy  disposition 
mechanism  becoming  increasingly  important  as  the  incident  energy  approaches  the  average 
thermal  energy  of  the  liquid  surface.  Lower  incident  energies,  of  course,  yield  proportionally 
lower  average  energy  transfer  just  because  the  energy  gap  between  the  initial  incident  energy  and 
the  average  thermal  energy  of  the  surface  is  smaller.  Nonetheless,  higher  rotational  excitations 
are  obtained  when  the  incident  energy  is  increased.  The  final  rotational  quantum  number 
distributions  corresponding  to  these  energy  distributions  are  shown  in  the  lower  panel  of  the 
figure  below.  Previous  experimental  work  by  Nesbitt  and  co-workers  yielded  comparable 
rotational  excitations  for  CO2  scattering  from  a  high  molecular  weight  molecular  liquid,  thus  our 
work  to  date  suggests  that  rotational  excitation  due  to  collisions  with  a  nonpolar  liquid  surface 
may  be  only  weakly  dependent  on  the  nature  of  the  liquid  species. 

This  project  is  ongoing,  with  an  investigation  of  the  dependence  of  the  rotational 
excitation  on  the  moment  of  inertia  of  the  incident  molecules  currently  in  progress,  as  is  a  study 
of  the  incident  angle  dependence  of  the  translation-to-rotation  energy  transfer.  Our  focus  here 
remains  on  a  full  characterization  of  the  collision  dynamics  in  a  prototypical  system  rather  than 
on  reproducing  the  results  of  a  particular  experiment  inasmuch  as  our  experience  to  date  has 
indicated  that  the  gross  features  of  energy  transfer  in  gas-liquid  collisions  are  interpretable  in 
terms  of  simple  kinematic  models. 


20 


%Occupaney 


40 


30 

20 

10 

0 


%Occupancy  vs.  %lnitial  Energy 


□  50kJ 

□  25kJ 


«  ~  V  (4)  a  , 

'  V  v  #  <%>  <£>  ^ 

%lnitial  Energy  N  ^ 


20 

> 

U  -  r 

c  15 

ns 


8 

o 


10 


%Occupancy  vs.  Rotational  Quantum  Number 


Rotational  Quantum  Number 


V  #  4?  4.  * 


Scientific  progress  and  accomplishments 

1.  For  both  atomic  and  molecular  liquids,  gas-surface  energy  transfer  is  dominated  by 
kinematic  effects.  Indeed,  the  gross  features  of  the  energy  transfer  obey  the  dependence 
on  masses  (in  particular,  the  ratio  of  the  impinging  species  mass  to  the  mass  of  the  liquid 
species)  proposed  by  Baule  to  describe  energy  transfer  at  a  solid  surface. 

2.  Energy  transfer  to  a  molecular  liquid  is  efficient  even  when  trapping  is  not. 

3.  In  the  case  of  liquid  nitromethane,  the  vibrational  degrees  of  freedom  of  the  nitromethane 
molecules  a  negligible  role  in  the  transfer  of  energy  due  to  collisions  (i.e.,  the  thermal 
accommodation  coefficient  differs  by  no  more  than  10%  when  vibrational  degrees  of 
freedom  are  “frozen”.) 

4.  For  a  molecular  impinging  gas,  collision  with  the  liquid  surface  leads  to  substantial 
rotational  excitation  of  the  species  reflected  from  the  surface.  Initial  results  suggest  that 
the  rotational  excitation  may  be  only  weakly  dependent  on  the  nature  of  the  liquid 
surface. 
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Development  of  Methods  for  Predicting  Solvation  and  Separation  for 
Energetic  Materials  in  Supercritical  Fluids 

Christopher  J.  Cramer,  Alek  V.  Marenich,  Adam  C.  Chamberlin,  Casey  P.  Kelly, 

Benjamin  J.  Lynch,  Ryan  M.  Olson,  Jason  D.  Thompson,  Paul  Winget,  and  Donald  G. 
Truhlar  (Department  of  Chemistry  and  Supercomputing  Institute,  University  of  Minnesota 
Minneapolis,  MN  55455) 

Statement  of  the  Problem  Studied 

Research  is  being  carried  out  at  the  University  of  Minnesota  on  two  projects.  First  is  the 
development  of  theoretical  methods  that  can  be  used  to  predict  the  solubilities  of  high-energy 
materials  in  various  fluids.  These  types  of  methods  have  a  variety  of  applications,  including  the 
development  of  improved  methods  for  the  safe  recycling/disposal  of  many  types  of  high-energy 
materials.  Second  is  the  development  of  new  methods  for  predicting  partial  atomic  charges  for 
modeling  the  crystallization  of  high-energy  materials  as  well  as  being  one  component  of  the 
solubility  calculations.  This  research  requires  the  development  of  charge  models  and  solvation 
models  that  may  be  used  for  a  wide  range  of  applications  including  the  modeling  of  crystalline 
materials.  During  the  MURI  grant  period,  we  have  developed  broadly  applicable  computational 
methods  for  solvation  free  energies  and  partial  atomic  charges.  Another  important  aspect  of  our 
research  has  been  improving  the  capability  and  portability  of  charge-model  and  solvation-model 
software  and  distributing  these  programs  (including  detailed  instruction  manuals  and  other 
related  documents)  to  fellow  MURI  researchers,  researchers  at  ARL,  and  to  the  scientific 
community  at  large.  We  are  implementing  them  into  easy-to-use  codes  that  are  freely  available 
via  the  World  Wide  Web.  Thus  Army  and  MURI  researchers  will  be  able  to  use  these  models 
and  methods  to  facilitate  the  design  of  practical  procedures  for  extraction,  recycling,  and  reusing 
materials. 

Since  the  MURI  project  began  in  2002,  we  have  focused  on  extending  existing 
continuum  solvation  models.  Eventually  it  would  be  desirable  to  have  a  model  for  supercritical 
solvents.  Our  work  so  far  has  been  focused  on  improving  the  performance  of  continuum 
solvation  models  for  predicting  solvation  free  energies  because  solubilities  can  be  computed 
using  solvation  free  energies  and  pure  solute  vapor  pressures.  Research  carried  out  under  this 
grant  examined  whether  calculated  solvation  free  energies  and  calculated  or  experimental  vapor 
pressures  could  be  used  to  make  accurate  predictions  of  solubility  (J.  D.  Thompson,  C.  J. 
Cramer,  and  D.  G.  Truhlar,  J.  Chem.  Phys.  119,  1661  (2003).  It  was  shown  that  in  most  cases, 
doing  this  leads  to  accurate  values  for  the  solubility,  even  for  solutes  that  are  solids  at  298  K. 
(The  vapor  pressure  can  be  calculated  using  our  continuum  solvation  models  because  it  is  related 
to  the  solvation  free  energy  of  the  solute  in  its  own  pure  liquid  phase,  or,  in  the  case  of  a  solid,  its 
supercooled  liquid  phase.  We  refer  to  this  free  energy  as  the  “self-solvation  free  energy.”)  Our 
continuum  solvation  models  (labeled  as  SMx  where  x  indicates  the  version  number)  separate  the 
solvation  free  energy  into  two  components:  the  long-range  bulk  electrostatic  contribution  arising 
from  a  self-consistent  reaction  field  treatment  using  the  generalized  Born  approximation  for 
electrostatics  is  augmented  by  the  non-electrostatic  contribution  arising  from  short-range 
interactions  between  the  solute  and  solvent  molecules  in  the  first  solvation  shell.  Our  continuum 
solvation  models  (i.e.,  solvation  models  for  ordinary  liquid  solvents)  utilize  the  generalized  Bom 
(GB)  model  in  a  self-consistent  reaction  field  step  to  account  for  long-range  electrostatic  effects. 
Note  that  the  GB  model  represents  the  solute  as  a  collection  of  atom-centered  spheres  and 
atom-centered  point  charges.  Our  model  also  accounts  for  solvation  phenomena  beyond 
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electrostatics,  including  cavitation,  dispersion,  hydrogen  bonding,  and  structural  changes  taking 
place  in  the  first  solvation  shell(s),  as  well  as  local  changes  in  the  permittivity  of  the  solvent,  by 
modeling  these  effects  as  being  proportional  to  the  solvent-accessible  surface  areas  (SASAs)  of 
the  atoms  in  the  solute.  The  constants  of  proportionality  are  called  atomic  surface  tension 
coefficients.  They  contain  parameters  (called  atomic  surface  tension  parameters)  that  are 
optimized  against  a  training  set  of  experimentally  known  free  energies  of  solvation  in  water  and 
organic  solvents  and  that  depend  on  a  set  of  experimental  solvent  descriptors  of  the  solvent. 
Progress  for  May  2002  -  August  2008  (Summary  of  the  Most  Important  Results) 

Improved  Charge  Models 

Throughout  the  MURI  project  we  have  been  continually  developing  and  improving  our 
charge  models,  charge  model  3  (CM3)  and  charge  model  4  (CM4),  that  provide  accurate  and 
reliable  charges  to  use  in  the  GB  model.  CM4  is  a  mapping  from  Class  II  charges  such  as 
Lowdin  (LPA)  or  redistributed  Lowdin  (RLPA)  for  diffuse  basis  sets  that  is  necessary  in  order  to 
correct  for  systematic  errors  in  electronic  population  analysis  techniques  for  charge-dependent 
observables.  Additionally  it  is  known  that  Lowdin  charges  are  not  stable  for  diffuse  basis  sets, 
and  because  diffuse  functions  may  be  important  for  an  even-handed  description  of 
conformational  energies  and  of  negative  functional  groups  found  in  many  of  the  compounds  of 
interest,  we  developed  a  new  population  analysis  method,  called  redistributed  Lowdin  population 
analysis  (RLPA).  While  Redistributed  LPA  charges  are  stable  for  diffuse  basis  sets,  they  still 
suffer  from  systematic  errors  similar  to  those  found  in  Lowdin  charges  for  non-diffuse  basis  sets. 
So  it  is  necessary  to  map  partial  atomic  charges  obtained  from  either  method  to  charges  which 
reproduce  experimentally  verifiable  properties  such  as  dipole  moments.  We  created  a  series  of 
improved  charge  models  ending  with  charge  model  4  (CM4).  CM4  uses  a  semiempirical 
charge-mapping  scheme  that  is  a  function  of  the  Mayer  bond  orders  in  the  molecule.  The  CM4 
charge  model  has  several  improvements  over  prior  class  IV  charge  models  developed  in  our 
research  group.  In  particular,  the  CM4  parameters  were  optimized  against  a  larger  and  more 
diverse  training  set  than  earlier  charge  models.  Second,  when  diffuse  basis  functions  are  used, 
CM4  maps  RLPA  charges  rather  than  Lowdin  population  analysis  charges.  In  addition,  the  CM4 
parameters  were  chosen  so  that  the  resulting  charges  are  not  as  sensitive  to  variations  observed  in 
the  Mayer  bond  order  when  diffuse  basis  functions  are  included  in  the  basis  set.  Lastly  the  C-H 
polarity  was  carefully  developed  using  a  specialized  set  of  19  different  compounds  designed  to 
constrain  the  C-H  bond  polarity  to  physical  reasonable  values.  An  additional  model  was  CM4M 
was  also  created.  While  both  CM3  and  CM4  were  designed  in  such  a  fashion  as  to  reliably 
reproduce  experimental  dipole  moments  independent  of  the  level  of  theory,  CM4M  was  designed 
to  best  reproduce  experiment  for  a  specific  class  of  DFT  functionals  namely  the  M06  series.  This 
allows  our  charge  model  to  take  full  advantage  of  a  very  reliable  and  accurate  level  of  theory. 
Validation  of  Charge  Models 

Our  development  of  the  charge  models  entailed  the  careful  validation  of  the  model  for 
predicting  accurate  charge  distributions  of  high-energy  materials  (HEMs)  and  compounds 
analogous  to  them.  The  charge  model  4  parameters  have  been  optimized  against  a  training  set  of 
dipole  moment  data,  and  validated  against  quadrapole  moments  computed  at  a  converged  level 
of  theory.  The  training  set  is  composed  of  416  compounds  containing  many  different  types  of 
functional  groups  that  one  encounters  in  organic  chemistry.  We  have  also  shown  that  atomic 
charges  calculated  with  CM4  yield  polarization  free  energies  (computed  from  the  GB  model)  for 
the  condensed  phase  that  are  less  sensitive  to  the  level  of  treatment  of  electron  correlation  than 
those  given  by  other  methods.  This  invariance  to  the  level  of  treatment  of  electron  correlation 
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demonstrated  by  the  CM4  model  has  been  critical  in  the  development  of  accurate  theoretical 
methods  for  modeling  the  solid-state  condensation  and  dissolution  of  HEMs.  CM4  is  an 
improvement  upon  CM3,  which  we  showed  accurately  reproduces  the  dipole  moments  of  HEMs. 
This  was  a  concern  because  the  functional  groups  found  in  HEMs  are  either  under-represented 
with  respect  to  other  functional  groups  present  in  the  training  set  (nitro  compounds)  or  not 
represented  at  all  (nitramines,  for  example)  in  the  CM3  training  set.  To  determine  whether  or  not 
the  CM3  parameters  are  applicable  for  predicting  accurate  charge  distributions  of  HEMs,  we 
assembled  a  test  set  of  compounds,  including  nitramide,  dimethylnitramine  (DMNA), 
1,3,3-trinitroazetidine  (TNAZ),  1 ,3,5-trinitro-v-triazine  (RDX),  and  hexanitrohexaaza- 
isowurtzitane  (HNIW),  and  compared  the  dipole  moments  for  different  conformations  of  these 
compounds  (14  in  all)  calculated  using  CM3  charges  to  high-level  theoretical  density  dipole 
moments.  Note  that  a  density  dipole  moment  is  calculated  from  the  one-electron  density  as  an 
expectation  value  of  the  dipole  moment  operator,  and  density  dipole  calculations  do  not  provide 
the  partial  atomic  charges  needed  for  condensed-phase  modeling.  These  tests  have  shown  that 
partial  atomic  charges  from  CM3  can  be  used  to  predict  dipole  moments  for  these  types  of 
compounds  with  similar  or  better  accuracy  as  for  compounds  in  the  CM3  training  set. 
Furthermore,  this  good  agreement  is  obtained  even  when  relatively  inexpensive  wave  functions 
are  used  for  the  solute,  which  is  important  for  larger  solutes  like  HNIW. 

Since  publishing  the  results  summarized  above,  we  have  been  in  close  contact  with  Betsy 
Rice,  who  has  begun  work  that  involves  testing  whether  CMx  (where  x  may  be  either  3  or  4) 
partial  atomic  charges  calculated  in  the  gas-phase  can  be  used  in  molecular  packing  calculations 
for  various  high-energy  materials.  Some  initial  calculations  that  use  CM3  gas-phase  partial 
atomic  charges  with  the  WMIN  lattice  minimization  program  have  yielded  encouraging  results. 
Improved  Solvation  Models  for  Room  Temperature 

The  accurate  reproduction  of  the  free  energy  of  solvation  depends  upon  two  factors, 
reproducing  the  electrostatic  properties  of  the  solute  and  reproducing  the  behavior  of  the  solvent. 
The  CMx  models  address  the  first  factor  however  the  latter  must  also  be  addressed.  To  this  end 
we  have  carefully  developed  a  series  of  solvation  models,  solvation  model  5.43R  (SM5.43R), 
solvation  model  6  (SM6),  and  solvation  model  8  (SM8),  wherein  we  have  explored  a  number 
means  to  improve  the  accuracy  of  the  model  with  the  ultimate  goal  of  making  our  model 
applicable  to  a  broad  range  of  solutes,  including  ions,  and  solvents.  In  particular  we  explored  the 
importance  of  including  explicit  water  molecules  near  highly  charged  sites  on  the  solute,  during 
the  development  of  solvation  model  6.  Additionally  we  explored  the  effect  of  the  Coulombic 
radii  upon  the  solvation  of  ions  in  nonaqueous  solvents  during  the  development  of  solvation 
model  8. 

The  SM8  model  improves  upon  earlier  SMx  universal  solvation  models  by  including  free 
energies  of  solvation  of  ions  in  nonaqueous  media  in  the  parametrization.  SM8  is  applicable  to 
any  charged  or  uncharged  solute  composed  of  H,  C,  N,  O,  F,  Si,  P,  S,  Cl,  and/or  Br  in  any 
solvent  or  liquid  medium  for  which  a  few  key  descriptors  are  known,  in  particular  dielectric 
constant,  refractive  index,  bulk  surface  tension,  and  acidity  and  basicity  parameters.  It  does  not 
require  the  user  to  assign  molecular-mechanics  types  to  an  atom  or  group;  all  parameters  are 
unique  and  continuous  functions  of  geometry.  It  may  be  used  with  any  level  of  electronic 
structure  theory  as  long  as  accurate  partial  charges  can  be  computed  for  that  level  of  theory;  we 
recommend  using  it  with  self-consistently  polarized  Charge  Model  4.  The  radii  used  for  aqueous 
solution  are  the  same  as  the  ones  developed  previously  for  the  SM6  aqueous  solvation  model, 
and  the  radii  for  nonaqueous  solution  vary  as  a  function  of  the  hydrogen  bond  acidity  of  the 


solvent.  The  non-electrostatic  terms  are  proportional  to  the  solvent-accessible  surface  areas  of 
the  atoms  of  the  solute. 

Validation  of  Solvation  Models  for  Room  Temperature 

These  models  were  then  carefully  optimized.  During  the  course  of  our  work  we  gradually 
improved  and  expanded  our  database  the  Coulomb  radii  in  organic  solvents  were  parameterized 
using  a  training  set  of  220  bare  ions  and  2 1  clustered  ions  in  acetonitrile,  methanol,  and  dimethyl 
sulfoxides,  and  the  radii  in  aqueous  solution  were  optimized  using  1 12  ionic  solutes,  and  31  ion- 
water  clusters.  The  non-bulk  electrostatic  terms  were  optimized  against  solvation  free  energies 
for  a  training  set  of  2346  solvation  free  energies  for  318  neutral  solutes  in  90  nonaqueous 
solvents  and  water  and  143  transfer  free  energies  for  93  neutral  solutes  between  water  and  15 
organic  solvents.  The  model  was  tested  with  three  density  functionals  and  with  four  basis  sets: 
6-31+G(d,p),  6-31+G(d),  6-31G(d),  and  MIDI16D.  The  SM8  model  achieves  mean  unsigned 
errors  of  0.5-0. 8  kcal/mol  in  the  solvation  free  energies  of  tested  neutrals  and  mean  unsigned 
errors  of  2.2-7. 0  kcal/mol  for  ions.  The  model  outperforms  the  earlier  SM5.43R  universal 
solvation  models  as  well  as  the  default  Polarizable  Continuum  Model  (PCM)  implemented  in 
Gaussian  98/03,  the  conductor-like  PCM  as  implemented  in  GAMESS,  Jaguar's,  continuum 
model  based  on  numerical  solution  of  the  Poisson  equation,  and  the  GCOSMO  model 
implemented  in  NWChem. 

Accounting  for  Temperature  Dependence 

Extraction  procedures  that  use  supercritical  solvents  are  carried  out  at  non-ambient 
temperatures  and  pressures.  In  fact,  the  ability  to  tune  the  properties  of  supercritical  solvents  by 
varying  the  temperature  and  pressure  is  one  of  their  key  advantages.  Thus  it  is  necessary  to 
extend  SM8  to  temperatures  other  than  298  K  in  order  to  carry  out  the  proposed  research. 
During  the  grant  period  we  developed  a  series  of  models  for  calculating  aqueous  solvation  free 
energies  as  a  function  of  temperature.  Solvation  model  6  with  temperature  dependence  (SM6T) 
used  the  SM6  for  the  computation  of  the  free  energy  of  solvation  at  room  temperature.  Since  the 
time  when  SM6T  was  created  we  developed  SM8.  The  aqueous  version  of  SM8  is  almost 
identical  to  the  aqueous  SM6  model,  and  it  has  the  same  Coulomb  radii  for  aqueous  solution. 
Thus  the  SM6T  parameters  work  equally  well  with  SM8  as  with  SM6.  Therefore  the  extension 
presented  here  is  made  using  SM8  rather  than  SM6.  We  refer  to  the  new  model  as  SM8T. 
Validation  of  Model  for  Temperature  Dependence 

An  important  step  that  is  necessary  in  order  to  develop  SM8T  is  the  creation  of  databases 
containing  accurate  experimental  data  for  a  diverse  set  of  molecules.  To  optimize  the  parameters 
in  SM6T,  we  used  a  database  of  aqueous  solvation  free  energies  for  solutes  containing  at  most 
the  elements  H,  C,  N,  O,  F,  S,  Cl,  and/or  Br  (a  total  of  4403  solvation  free  energy  data  points  for 
348  solutes).  Because  the  SMxT  (x  =  6  or  8)  models  are  the  first  theoretical  models  of  their  type 
for  predicting  temperature-dependent  solvation  free  energies,  the  development  of  these  models 
was  an  important  step  in  determining  the  accuracy  to  which  temperature-dependent  solvation 
free  energies  can  be  predicted.  When  tested  against  the  database  described  above,  which  includes 
aqueous  solvation  free  energies  in  the  temperature  range  of  273  K  to  373  K  (corresponding  to  the 
freezing  and  boiling  points  of  pure  water),  SM8T  achieves  a  root  mean  squared  error  of  0.09 
kcal/mol.  If  one  were  to  assume  that  temperature  has  no  effect  on  the  aqueous  solvation  free 
energy  (null  hypothesis),  the  root  mean  squared  error  would  be  significantly  higher  (0.70 
kcal/mol).  The  partial  atomic  charge  models  developed  for  and  used  in  the  solvation  models  can 
also  be  used  in  studies  of  the  crystal  structure  and  energetics  of  high-energy  materials. 
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Method  Development  and  Simulations  of  Physical  Processes  and  Chemical 

Reactions  in  the  Various  Phases 

Donald  L.  Thompson  (University  of  Missouri  -  Columbia) 

The  focus  of  the  research  at  the  University  of  Missouri  was  on  practical  methods  for 
simulations  and  rate  calculations  of  physical  and  chemical  processes  in  the  gas,  liquid,  and  solid 
phases  of  energetic  materials;  more  specifically: 

•  Development  of  accurate  methods  for  simulating  melting  and  predicting  melting  points  of 
organic  solids. 

•  Development  practical  methods  for  computing  reaction  rates  of  reactions  in  condensed  phases. 

•  Development  of  methods  for  fitting  ab  initio  potential  energy  surfaces  and  for  performing 
direct  dynamics  simulations. 

•  Using  quantum  chemistry  methods  to  determine  the  reaction  pathways  for  unimolecular  and 
|  bimolecular  reactions  of  energetic  molecules. 

Melting  of  Energetic  Materials 

A  significant  result  of  our  work  is  the  development  and  demonstration  of  methods  for 
simulating  solid-to-liquid  phase  transitions  in  energetic  materials  to  predict  their  melting  points. 
Prior  to  this  grant,  working  with  Drs.  Betsy  M.  Rice  (ARL)  and  Dan  C.  Sorescu  (Okla.  St.  U.), 
we  developed  crystal  models  for  organic  materials  that  accurately  describe  a  wide  range  of 
energetic  materials.  In  the  current  grant  we  have  used  these  models  in  studies  of  melting 
transitions  of  energetic  materials.  The  simulation  of  melting  and  the  prediction  of  melting  points 
have  been  considered  difficult  problems.  The  difficulty  stems  from  the  free  energy  barrier  for 
the  formation  of  a  liquid-solid  interface;  which  can  cause  superheating  in  a  perfect  crystal  and 
thus  an  overestimation  of  the  melting  point  of  the  real  material.  We  have  investigated  various 
ways  of  performing  the  simulations  to  determine  the  most  practical  one  for  use  in  studies  of 
energetic  materials.  Prior  to  our  work,  the  simulations  had  been  restricted  to  atomic  solids  such 
as  rare  gases  and  metals.  We  first  carried  out  a  series  of  calculations  for  argon  to  test  methods, 
and  then  we  extended  the  studies  to  energetic  crystals. 

The  introduction  of  vacancies  in  a  crystal  eliminates  the  free  energy  barrier  that  causes 
superheating  in  a  perfect  crystal.  Alternatively,  the  simulations  can  be  initiated  with  the  solid  in 
contact  with  liquid.  Molecular  dynamics  studies  of  melting  of  nitromethane  using  the 
intermolecular  interaction  potential  of  Sorescu  et  al.  (SRT)  [D.  C.  Sorescu,  B.  M.  Rice,  and  D.  L. 
Thompson,  J.  Phys.  Chem.  B  104,  8406  (2000)]  have  been  carried  out  by  two  methods:  (1)  void- 
nucleated  melting  with  the  gradual  heating  of  the  lattice  and  (2)  equilibration  of  coexisting  liquid 
and  solid  phases.  The  results  are  in  near  agreement  with  each  other;  the  small  difference  is 
attributed  to  the  hysteresis  effect  associated  with  the  direct  heating  process.  The  computed 
values  of  the  melting  temperature  are  in  good  agreement  with  the  experimental  data  at  various 
values  of  pressure  ranging  from  1  atm  to  30  kbar.  Since  the  void-induced  melting  approach  is 
much  easier  to  use,  we  have  adopted  it  as  the  practical  approach  for  predicting  melting  points  of 
energetic  materials.  Following  this  work,  which  used  our  well-tested  force  field  for 
nitromethane,  we  extended  the  studies  to  practical  energetic  materials,  which  provide  critical 
tests  of  the  models  and  simulation  methods.  We  find  that  the  results  are  sensitive  to  details  of 
the  potential  and  that  for  large,  complex  molecules  the  simulations  are  somewhat  more 
demanding;  however,  the  results  show  that  approaches  that  we  are  using  will  be  effective  and 
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practical  for  predicting  melting  points.  This  work  was  done  in  collaboration  with  Dr.  Betsy  Rice 
at  ARL. 

We  have  developed  methods  for  simulations  of  the  melting  transitions  in  perfect  and 
imperfect  crystals  of  energetic  solids.  The  focus  was  on  developing  and  testing  methods,  and 
thus  most  of  the  simulations  were  performed  for  nitromethane  because  of  its  relatively  simple 
molecular  structure  and  the  availability  of  the  accurate  SRT  force  field;[D.  C.  Sorescu,  B.  M. 
Rice,  &  D.  L.  Thompson,  J.  Phys.  Chem.  B  104,  8406  (2000).]  however,  studies  of 
dimethylnitramine  (DMNA),  trinitroazetidine  (TNAZ),  and  hexahydro-1  ,3,5-trinitro- 1 ,3,5-s- 
triazine  (RDX)  were  carried  out.  We  have  shown  that  these  methods  can  be  used  to  accurately 
predict  the  thermodynamic  melting  point  and  superheating  temperatures,  and  we  have  used  them 
to  investigate  the  molecular-level  details  of  the  melting  transitions  in  these  materials.  This  is 
important  in  the  development  of  realistic  models  of  the  chemistry  energetic  materials  because 
phase  transitions,  particularly  melting,  are  precursory  to  chemical  initiation.  Thus,  our  studies 
not  only  validate  the  force  field  models  but  also  provide  physical  insights  into  the  changes  that 
occur  when  energetic  materials  are  subjected  to  heat  or  pressure.  An  illustration  of  the  accuracy 
of  the  methods  is  given  in  the  following  figure,  which  shows  the  computed  and  measured 
melting  points  of  nitromethane  as  functions  of  pressure. 


The  blue  circles  show  the  results  for  melting  initiated  at  voids  within  the  crystal  [P.  M.  Agrawal,  B.  M. 
Rice,  &  D.  L.  Thompson,  J.  Chem.  Phys.  119,  9617  (2003).]  and  the  red  triangles  are  the  computed 
melting  points  for  the  crystal  in  contact  with  Ar  gas  [A.  Siavosh-Haghighi  and  D.  L.  Thompson, 
manuscript  in  preparation.].  The  experimental  results  [G.  J.  Piermarini,  S.  Block,  &  P.  J.  Miller,  J.  Phys. 
Chem.  93,  457  (1989)  &  W.  M.  Jones  &  W.  F.  Giauque,  J.  Am.  Chem.  Soc.  69,  983  (1947).]  are 
represented  by  the  green  diamonds;  the  solid  curve  is  a  least-squares  fit  to  experimental  values. 
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Practical  Rate  Calculation  Methods 

While  the  main  thrust  of  the  work  has  been  explicit  MD  simulations  of  processes  in 
condensed  phase  energetic  materials,  we  also  began  development  of  methods  for  accurate  rate 
calculations  that  complement  the  results  of  simulations,  especially  methods  that  can  be  used  to 
either  save  computer  time  or  to  treat  reactions  that  cannot  be  dealt  with  in  a  simulation. 
Specifically,  we  focused  on  methods  for  computing  rates  in  the  liquid  phase  for  the  initial  efforts 
to  arrive  at  methods  for  all  condensed  phases.  In  our  initial  effort,  we  used  relatively  simply 
reactions  and  performed  the  studies  for  reactions  in  simple  solvents;  however,  the  approach  can 
be  extended  to  more  complex  systems  of  practical  interests  in  energetic  materials  research. 

The  solvent  effects  on  rates  of  chemical  reactions  have  been  extensively  studied  over  the 
years.  Much  of  the  theoretical  work  is  based  on  stochastic  methods  such  as  the  Langevin  theory. 
These  methods  allow  for  computations  of  rate  constants  without  explicit  consideration  of  the 
detailed  dynamics  of  the  solvent,  and  are  thus  powerful  for  treating  complex  systems  where  full¬ 
dimensional  dynamical  calculations  are  too  costly.  For  instance,  in  the  Langevin  theory,  the 
solvent  influence  on  the  reacting  particle  is  modeled  by  a  random  fluctuating  force  and  a  friction 
kernel.  Although  these  stochastic  methods  are  useful  for  qualitatively  understanding  reactions  in 
condensed  phases,  rate  constants  cannot  be  directly  obtained  for  a  given  realistic  system  because 
quantities  involved  in  the  models  concerning  the  solvent  effects  are  generally  unknown  and  must 
be  determined  by  examining  the  detailed  dynamics  of  the  system.  Molecular  dynamics  (MD) 
simulations  provide  a  useful  means  for  realistically  modeling  detailed  dynamics  and  elucidating 
reaction  mechanisms.  However,  a  full-dimensional  MD  study  is  computationally  expensive  and 
may  not  be  feasible  for  complex  systems.  An  alternative  approach  is  to  map  the  real  system  to  a 
stochastic  model  by  performing  a  few  MD  calculations  to  determine  the  fitting  parameters  in  the 
model.  We  have  explored  such  an  approach. 

The  basic  idea  is  that  the  reacting  molecule  is  subjected  to  random  collisions  with  a  mean 
frequency  a.  Between  collisions,  the  dynamics  of  the  molecule  is  governed  solely  by  the 
equations  of  motion  for  the  molecule  itself.  The  effects  of  the  solvent  are  modeled  by  random 
collisions.  Under  the  strong  collision  assumption,  the  velocities  of  the  molecule  are  randomized 
after  each  collision  by  sampling  from  the  Boltzmann  distribution  while  the  coordinates  are  held 
fixed.  Since  only  the  motion  of  the  molecule  is  explicitly  treated,  the  method  is  far  less  costly 
than  full-dimensional  MD  simulations.  The  parameter  used  to  model  the  solvent  effects  is  the 
collision  frequency  a,  which  is  related  to  the  liquid  density  p  of  the  system.  There  have  been 
studies  where  both  MD  simulations  and  the  stochastic  dynamics  method  were  employed  on  the 
same  system;  but  no  attempt,  to  our  knowledge,  had  previously  been  made  to  quantitatively 
relate  a  and  p  for  any  system. 

We  first  investigated  this  idea  with  an  application  to  cis-trans  isomerization  in  HONO, 
which  is  one  of  the  simplest  molecules  that  exhibit  cis-trans  isomerization  and  one  that  has  been 
extensively  studied  experimentally  and  theoretically. [Y.  Guo  and  D.  L.  Thompson,  J.  Chem. 
Phys.  120,  898  (2004).]  We  used  a  model  system  with  a  HONO  molecule  embedded  in  liquid 
krypton  and  performed  MD  simulations  to  obtain  isomerization  rates  at  several  liquid  densities. 
We  also  computed  rate  constants  using  the  stochastic  dynamics  method  for  a  wide  range  of 
collision  frequencies  a  for  these  two  sets  of  computed  rates.  We  then  examine  the  relationship 
between  the  liquid  density  p  and  the  collision  frequency  a.  Comparisons  of  the  two  sets  of  the 
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computed  rates  show  that  for  a  wide  range  of  liquid  densities  there  is  a  simple  linear  relation 
between  the  liquid  density  p  and  the  collision  frequency  a,  that  is,  a  =  cp.  This  suggests  that 
once  the  constant  c  is  determined  from  a  molecular  dynamics  calculation  at  a  single  density,  the 
reaction  rates  can  be  obtained  from  stochastic  dynamics  calculations  for  the  entire  range  of  liquid 
densities  where  a  =  cp  holds.  The  applicability  of  the  combined  molecular  dynamics  and 
stochastic  dynamics  approach  provides  a  practical  means  for  obtaining  rate  constants  at 
considerable  savings  of  computer  time  compared  to  that  required  by  using  full-dimensional 
molecular  dynamics  simulations  alone.  Given  the  success  of  this  study,  we  next  carried  out  a 
similar  study  of  bond  fission  in  HONO  (i.e.,  HONO  — »  OH  +  NO).[Y.  Guo  and  D.  L.  Thompson, 
Chem.  Phys.  Letters  418,  351-354  (2006).] 

Intuitively,  the  solvent  effects  on  dissociation  reactions  should  be  different  from 
those  for  isomerization;  the  cage  effect  due  to  the  solvent  that  forces  recombination  should  be  a 
prominent  factor  in  dissociation  reactions.  Thus,  the  relation  between  a  and  p  as  well  as  the 
dependence  of  the  reaction  rate  on  p  are  expected  to  be  different  from  that  for  isomerization.  To 

establish  the  relation  between  a  and  p,  we  have  performed  two  sets  of  calculations.  We 
calculated  the  rate  as  a  function  of  a  using  the  stochastic  dynamics  method.  We  also  calculated 
the  rate  at  four  liquid  densities  by  using  the  full-dimensional  MD  simulations  with  a  HONO 
molecule  inside  liquid  Kr.  The  relation  between  a  and  p  was  found  by  pairing  a  and  p  that 
have  the  same  rate.  In  our  study  of  the  cis-trans  isomerization  of  HONO  in  liquid  Kr,  we  found 
that  a  and  p  obey  a  simple  linear  relation,  a  =  cp  .  For  the  O-N  bond  dissociation  of  HONO  in 
liquid  Kr,  we  found  that  the  relationship  is  more  complicated,  but  is  accurately  described  by  an 
analytical  expression.  We  have  also  found  that  the  solvent  effects  are  different  for  the 
isomerization  and  dissociation  reactions.  For  isomerization  the  rate  is  in  the  low-collision  regime 
and  increases  with  increasing  liquid  density,  whereas  for  dissociation  the  rate  is  in  the  high- 
collision  regime  and  decreases  with  increasing  liquid  density. 

In  both  cases,  based  on  one  MD  simulation  at  a  single  density,  the  reaction  rates  can  be 
obtained  from  stochastic  dynamics  for  a  wide  range  of  liquid  densities.  The  approach  thus 
provides  tremendous  savings  of  computational  time  compared  to  full-dimensional  MD 
simulations.  Future  work  should  be  carried  out  to  extend  the  studies  to  larger,  more  complicated 
systems,  including  more  realistic  energetic  molecules  and  in  reactions  in  neat  liquids. 

Reaction  Pathways  for  Unimolecular  Dissociations  of  Energetic  Molecules 

A  major  challenge  to  theoretical  chemists  and  a  particular  need  in  energetic  materials 
research  is  the  development  of  methods  for  simulations  and  rate  calculations  for  sequential, 
branching  reactions.  The  practical  problems  range  from  soot  formation  to  the  decomposition  of 
an  energetic  material  such  as  RDX.  We  have  pursued  two  lines  of  research  on  this  problem. 
First,  we  used  quantum  chemistry  to  explore  the  transition  states  and  intermediates  in  the 
chemical  decomposition  of  energetic  molecules.  The  other  focused  on  methods  for  fitting  ab 
initio  potential  energy  surfaces  and  for  performing  direct  dynamics  simulations.  Quantum 
chemistry  studies  of  the  decomposition  pathways  of  1,3,3-trinitroazetidine  (TNAZ)  [S.  Alavi,  L. 
M.  Reilly,  and  D.  L.  Thompson,  J.  Chem.  Phys.  119,  8297  (2003).]  and  dimethylnitramine  [G.  F. 
Velardez,  S.  Alavi,  and  D.  L.  Thompson,  J.  Chem.  Phys.  123,  074313  (2005).]  have  been  carried 
out.  Possible  decomposition  transition-states,  intermediates,  and  products  of  TNAZ  were 
identified  and  the  structures,  energies,  and  vibrational  frequencies  were  determined  at  the 
B3LYP/6-3  1G(J,/i)  level.  Four  major  pathways  are  apparent.  Two  pathways  are  initiated  by  the 


29 


fission  of  the  N-NO2  and  C-NO2  bonds  to  yield  radical  intermediates,  while  the  other  two 
pathways  involve  the  molecular  elimination  of  HONO.  Energy  profiles  for  the  pathways  and 
possible  routes  to  some  of  the  experimentally  observed  species  of  TNAZ  decomposition  have 
been  determined.  The  energies  required  to  initiate  the  NO2  bond  fission  pathways  are  4  to  8 
kcal/mol  lower  than  the  HONO  elimination  pathways.  In  the  gas  phase,  the  NO2  elimination 
pathways  will  be  the  dominant  routes  for  TNAZ  decomposition.  In  the  condensed  phase, 
however,  this  trend  may  be  reversed. 

Examination  of  the  reported  experimental  and  theoretical  results  from  several  researchers 
shows  that  there  is  an  significant  dispersion  in  the  values  of  the  Arrhenius  parameters  of  the  rate 
constants  for  the  unimolecular  decomposition  of  DMNA.  The  structures  and  energies  of  the 
reactant,  intermediates,  products,  and  transition  states  of  the  initial  stages  of  DMNA 
decomposition  of  have  been  determined  by  quantum  chemical  calculations  at  several  levels  of 
theory.  Kinetic  parameters  for  these  decomposition  steps  have  been  obtained  using  the  RRKM 
formalism  for  the  range  300-5000  K.  The  pathways  considered  are  NO2  elimination,  HONO 
elimination,  and  NNO2-NONO  rearrangement.  The  NO2  elimination  channel  is  the  main  channel 
of  gas  phase  decomposition  of  DMNA  in  the  300-5000  K  range  of  temperatures.  The  HONO 
elimination  channel  has  the  next  lowest  decomposition  energy  but  at  high  temperatures  the  nitro- 
nitrite  rearrangement  competes  with  it. 

The  rate  constants  as  functions  of  temperature  for  NO2  elimination,  HONO  elimination, 
and  nitro-nitrite  arrangement  have  been  obtained  from  values  of  the  energies  and  frequencies 
using  RRKM  theory;  Arrhenius  plots  are  shown  in  the  figure  below,  where  the  curves  are 
identified  as  follows:  ( — )  NO2  elimination,  ( — )  HONO  elimination,  and  (•••)  nitro-nitrite 
rearrangement. 


l/T  (K'1) 


The  Arrhenius  parameters  for  NO2  elimination  are  in  good  agreement  with  the  values  of 
previous  theoretical  studies.  The  HONO  elimination  channel  has  a  lower  value  of  the  pre¬ 
exponential  factor  compared  to  NO2  elimination  because  the  transition  state  for  the  HONO 
elimination  pathway  is  entropically  unfavorable.  Clearly,  the  NO2  elimination  is  the  main 
decomposition  channel  of  DMNA  in  the  range  of  300-5000  K.  The  HONO  elimination  has  a 
higher  rate  than  the  nitro-to-nitrite  rearrangement  but  at  high  temperatures  the  rates  of  these  two 
reactions  becomes  comparable.  The  experimental  evidence  is  that  the  presence  of  NO  in  the 


30 


DMNA  decomposition  products  is  due  to  the  nitro-nitrite  rearrangement  followed  by  breaking  of 
the  weak  NO-NO  bond.  Our  results  show  that  at  low  temperatures,  this  channel  is  the  third  in 
importance  and  becomes  comparable  with  the  HONO  elimination  channel  only  at  higher 
temperatures. 

We  carried  out  quantum  chemistry  computations  to  investigate  the  reactions  of  the  radicals 
produced  in  the  decomposition  of  cyclic  nitramines. [Methylene  amidogen  H2CN  is  an 
intermediate  in  the  thermal  decomposition  of  HDX  and  RDX  and  in  a  number  of  other  chemical 
processes.[G.  F.  Adams  and  R.  W.  Shaw,  Jr.  Annu.  Rev.  Phys.  Chem.  1992,  43,  311.]  Recently, 
Nizamov  and  Dagdigian  [B.  Nizamov  and  P.  J.  Dagdigian,  J.  Phys.  Chem.  A  2003,  107,  2256.] 
studied  room-temperature  chemical  reactions  involving  H2CN  and  reported  the  room- 
temperature  rate  constant  for  the  H2CN  +  OH  reaction.  Their  results  indicate  that  the  hydrogen 
abstraction  channel  to  form  HCN  +  H20  is  the  predominant  reaction  channel.  We  performed 
high-level  ab  initio  calculations  to  identify  the  likely  products  of  the  ground  state  reaction 
between  H2CN  and  OH,  with  the  eventual  goal  of  evaluating  rate  constants  for  these  reactions. 
The  reaction  channels  considered  are: 


H2CN  +  OH  -> 

HCN  +  H20 

(1) 

-> 

h2cnoh 

(2) 

-> 

h2conh 

(3) 

-> 

H2CN(H)0 

(4) 

-> 

ch3no 

(5) 

-> 

CH30N 

(6) 

-> 

HCON  +  H2 

(7) 

-> 

HCNO  +  H2 

(8) 

-» 

HNCO  +  H2 

(9) 

-> 

h2cnh  +  0 

(10) 

All  possible  isomeric  forms  of  each  product  were  investigated.  With  the  exception  of  HCN, 
H20,  and  H2,  both  singlet  and  triplet  states  were  examined  for  the  molecular  products.  We  have 
|  determined  transition  states  and  for  transformations  between  the  possible  products. 

More  Efficient  use  of  Ab  Initio  Potential  Energy  Surfaces 

An  important  part  of  predicting  the  behavior  of  energetic  materials  is  better  practical 
methods  for  computing  the  rates  of  chemical  reactions.  A  critical  problem  is  the  formulation  of 
PESs.  For  purely  predictive  methods  it  necessary  that  the  PESs  be  based  on  quantum  chemistry 
calculations,  but  this  remains  a  major  obstacle  to  applications  to  large  polyatomic  molecules  and 
radicals.  The  most  straightforward  way  to  make  use  of  ab  initio  results  is  to  simply  compute  the 
forces  “on  the  fly”  during  a  simulation;  however,  the  difficulty  here  is  that  studies  are  restricted 
to  the  lower-level  quantum  methods,  which  are  often  not  sufficiently  accurate  for  reliable 
predictions  of  rates.  The  common  alternative  to  “direct  dynamics”  is  to  fit  the  ab  initio  energies 
to  an  analytical  function,  which  allows  for  much  more  rapid  evaluations  of  the  forces  during 
simulations  as  well  as  scaling  of  the  ab  initio  results  to  correct  for  the  inaccuracies  inherent  in  a 
low-level  theory.  Fitting  a  global  PES  is  extremely  tedious  and  not  readily  generalizable,  thus  a 
huge  investment  of  human  labor  is  required  in  each  case.  Furthermore,  the  analytical  forms  are 
not  sufficiently  flexible  to  accurately  fit  the  ab  initio  points.  This  was  our  motivation  for 
suggesting,  in  the  1970’s,  that  local  fitting  functions  should  be  used  to  fit  ab  initio  energies.  We 
proposed  using  cubic  splines,  which  provide  the  desired  flexibility  to  yield  accurate,  smooth  first 
and  second  derivatives  but  requires  a  fairly  dense  grid  of  ab  initio  points.  [  D.  R.  McLaughlin  and 
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D.  L.  Thompson,  J.  Chem.  Phys.  59,  4393  (1973).  The  modified  Shepard  methods  suggested  by 
Ischtwan  and  Collins  [J.  Ischtwan  and  M.  A.  Collins,  J.  Chem.  Phys.  100,  8080  (1994).]  in  1994 
are  superior  to  cubic  splines.  They  have  combined  surface  fitting  with  MD  simulations  in  an 
iterative  scheme  for  successively  and  automatically  improving  the  PES.  This  procedure  selects 
the  locations  for  additional  ab  initio  calculations  in  regions  of  configuration  space  that  are 
dynamically  important.  The  procedure  is  simple  and  readily  automated  but  gradients  and 
Hessians,  which  are  not  readily  available  in  the  high-level  ab  initio,  are  required  at  every  point. 
We  have  developed  interpolative  moving  least-squares  (IMLS)  fitting  procedures  that  promise  to 
be  more  useful  than  other  methods,  mainly  because  it  does  not  require  derivatives. [Y.  Guo,  L.  B. 
Harding,  A.  F.  Wagner,  M.  Minkoff,  and  D.  L.  Thompson,  J.  Chem.  Phys.  126,  104105  (2007).] 

We  have  begun  applications  of  the  IMLS  methods  to  reactions  important  in  the 
decomposition  of  nitramines,  specifically  a  study  of  the  unimolecular  decomposition  of  H2CN, 
which  is  currently  under  study  in  Professor  Paul  Dagdigian’s  laboratory  at  Johns  Hopkins.  Ab 
initio  energies  were  fit  and  rates  of  dissociation  computed.  These  results  are  available  in 
advance  of  Professor  Dagdigian’s  measurements  of  the  rates,  thus  they  will  be  truly  predictive. 

The  first  objective  was  to  establish  the  accuracy  of  this  approach  for  a  significant  chemical 
reaction.  The  long-term  goal,  which  is  beyond  the  scope  of  this  project,  is  to  develop  software 
that  manages  the  quantum  chemistry  calculations  in  direct  response  to  the  needs  of  a  MD 
simulation  to  provide  an  ab  initio  prediction  of  a  specified  reaction  rate.  If  successful,  and  the 
prospects  are  good,  this  could  make  MD  studies  of  reactions  a  practical  tool  for  all  who  are 
interested  in  reaction  kinetics  rather  than  a  technique  limited  to  use  by  experts. 
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